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ABSTRACT 


The Finite Element Analysis Program (FEAP) was expanded 
to solve linear and nonlinear, two and three dimensional 
heat conduction problems. The usual types of , boundary 
conditions, including radiation, may be specified. A wide 
range of two- and three-time level schemes for the solution 
of time dependent problems is available In the program and a 
discussion of those most commonly used is presented. The 
algorithms for the solution of typical practical problems are 
described and several numerical examples are presented. The 
results are compared with the available analytical 
solutions. A listing of this expanded version of FEAP and 


the corresponding user's instructions are provided. 
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I. INTRODUCTION 


This thesis describes an application of the Finite 
Element method to Heat Transfer analysis through the use of 
a computer program. This program was designed to solve 
linear and nonlinear, steady and unsteady, two and three 
Miensional heat conduction problems’ involving temperature 
dependent thermophysical properties and complicated 
radiation/convection boundary conditions. 

The Finite Element Analysis Program (FEAP), programmed 
by Professor R. L. Taylor in the Department of Civil 
Engineering of the University of California, Berkeley, 
served as the point Sr pantan for the present code. It 
has now been expanded with two additional modules for Heat 
Transfer analysis and time integration of first order 
equations, respectively, but the original characteristics of 
a research and educational tool in which the various modules 
can be changed or added to as desired, were maintained. 

Although all equations are retained in core, the program 
can handle realistic engineering problems with several 
hundred unknowns in most computer systens. 

While the algorithms to solve steady heat transfer 
problems are well studled and defined, more research has to 
be done on the solution of unsteady problems. The 
Zienkiewicz two- and three-time level schemes were 


integrated in the program in such a way that a wide range of 





choices of time integration algorithms is available to the 
user. A limited study was performed in order to determine 
the characteristics of the most used members of those 
families of numerical schemes. 

The FEAP code is written in FORTRAN IV language and this 
version was constructed and tested on an IBM 360/67 computer 
system with 0S/67 release 18. Other systems and 


installations wlll require modifications. 





Il. STATEMENT OF OBJECTIVES 


A. EQUATION OF CONDUCTION. INITIAL AND BOUNDARY 
CONDITIONS 
The problem considered for solution is thoroughly 
developed by Arpaci in Ref. 1 and is mathematicaly 
described in the region 2 by the equation 
oT 


Pac See VE VID cQ (1) 


subjected to boundary conditions 


T = T on Γι (2) 
W 
and 
πα. E uus 0 on > (3) 
and initial condition 
lim T = T (4) 
t>0 0 
V is the gradient operator, I| and T2; are mutually 


exclusive parts of the boundary of the region 2 , T is the 
temperature and t is the time. The thermal Bay pc, the 
thermal conductivity K and the rate of internal heat 
generation Q are thermophysical properties dependent on 
temperature, 

In equation (3) A is an outward unit vector normal to 


the boundary surface whlle q, q represent respectively 


c^ Ar 
the imposed heat flux and the rates of heat flow per unlit 


area due to convection and radiation defined as 


q_ = hi (T - T. ) | (5) 


C C 


and 





4 


4 
qp = Fo (T - l 


) = h. ( T - p ) (6) 

In (5) hc is the temperature dependent convective heat 
transfer coefficient. in (6) the parameter h. is defined by 
the expression 

Е 2 2 Е 
N το ( T7 + er ΠΠ d ), where F is 

the radiant exchange factor апа со the Stefan-Boltzman 
constant. Tac and nr are the equilibrium temperatures for 


which, respectively, no convection and radiation occurs. 


B. NUMERICAL SOLUTION, DISCRETIZATION BY THE GALERKIN 
METHOD. MATRIX FORMULATI ON 
The spacewlse discretization of .equation (1) using 
cartesian coordinates can be acomplished by Galerkin's 
principle as shown by Zienkiewicz In Ref. 2 and Lew in Ref. 
B. 
Let the unknown function T be approximated, throughout 
the solution domaln at any time t, by the relationship 
T - EN Gy) T; (t) = <N> {T} (7) 
where N, are the usual shape functions defined piecewise 
element by element, T: being the nodal parameters. 
The result is 
[K] {T} + [C] {T} + {F} = (0) (8) 
where [K] and [C] are symmetric matrices defined, on the 


element level, as 


T T 
e _ ÓN ÓN ÓN ON 
[K]~ = Јое (<> eK “су πι ky + 
«SN; ¿NS k ) dQ + 
Š Z Š Z 7,2 
Ге ( h. + h. ) <N> <N> dT (9) 
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[C] ° = NC <N> de (10) 


The load vector F is 
{F}® = -f _<nd Q da + 
f ρον (q - hT,, - h,T,, ) dr (11) 

where the surface integrals are performed only over the 
surfaces where the prescribed boundary condition applies. 

lt must be noted that the set of equations (8) is 
nonlinear since the matrices [K] , [C] and vector {Flare 
dependent on T. 

The system of equations (8) may now be solved by any 
automated numerical technique. Code FEAP was written for 


this purpose and is described in the rest of this thesis. 
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iii. FINITE ELEMENT ANALYSIS PROGRAM 





The Finite Element Analysis Program (FEAP) was 
programmed and published by Taylor in Ref. 4& where it is 
well explained. The reader should consult this reference for 
details about it. 

This thesis is an expansion of the original program 
towards the solution of Heat Conduction problems and as such 
only the main features of the basic program will be referred 
in this section. A more detailed explanation is reserved for 
the new subroutines added to the program. The complete 
user's instructions for FEAP are given in Appendix A and 
they complement this section. 

A. MAIN PROGRAM 

FEAP can be separated into two basic parts: 

a) Data input module and preprocessor 

b) Solution and output modules 

The data input module must transmit sufficient 
information to the other modules so that each problem can be 
solved. All the input data is stored in a single array 
(integer array M) which is partitioned to store all the data 
arrays, as well as the global arrays, e.g., stiffness, mass, 
load, etc, 

The total capacity of the program is controlled by the 


dimension of the array M in the blank common of the main 
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program and the corresponding value of the variable MAX. 

The partition of M is performed in the control routine 
(subroutine PCONTR) and the data used for it is supplied by 
the user in the title and control information cards (first 
two cards for every run of the program). 

Then the nodal coordinates, element connections, 
material properties, nodal loading and boundary conditions 
may be input using the macro control statements. They are 
interpreted by subroutine PMESH and each control is a 
function independent of the others. 

To clarify what was said till this point, the discussion 
of a simple example may be useful. Consider the mesh 
presented in Figure land assume that quadratic rectangular 
elements as the one defined in Figure 2 are used. This 
example is completly solved in Appendix B and will be used 
throughout this text. 

The sector of a circular ring was divided in six 
elements and the 33 nodes were numbered. The order of 
numbering is not crucial but in order to improve the profile 
of non-zero coefficients the following general rule should 
be used. 

The numbering should be such as to minimize the nodal 
difference for each element (maximum node nummber minus 
minimum node number). 

The user can now proceed to the preparation of data for 
the program. The first step consists of specifying the 


problem title and the control information. The latter is 


¡Figures are grouped at the end pp 47-79 
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clearly 33 nodes, six elements, one material set, two 
spatial dimensions and eight nodes per element. !f this is a 
Heat Transfer problem we have one unknown per node 
(temperature). 

The program expects now the data cards for mesh 
description, Any analysis require at least 

a) coordinate data which follows the macro command COOR 

b) connectivity table which follows the macro command 
ELEM 

c) material data which follows the macro command MATE 

In addition, most analysis will require specification of 
nodal boundary restraint conditions, macro BOUN, and the 
corresponding nodal force or displacement values, macro 
FORC. The term displacement refers to the value of the 
unknown for each specific problem. іп a Heat Transfer 
problem the unknown Is temperature while in an Elasticity 
analysis it Is a geometrical displacement. 

The TEMP macro command (temperature data) shall not be 
used in a Heat Transfer analysis to specify nodal 
temperature. It may be used to specify auxilary nodal 
quantities. 

If in our problem we assume the line defined by nodes 
31, 32 and 33 at the temperature of 20, the line defined by 
nodes 1, 2 and 3 at 820 and Insulated elsewhere, the macro 
control statements and correspondent data cards will be 


COOR 
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31 . É 
+ 5.116666667 


29 «2833332232 
2 э 

52 .3 

3 5 œl 

35 <2 


5 5.116666667 
30 283333333 
(blank card) 
POLA 
1 33 1 
(blank card) 
ELEM 
1 1 1 6 
(blank card) 


BOUN 
1 1 =1 
5 = 
31 1 sI 
33 al 


(blank card) 
FORC 





(blank card) 
MATE 

1 2 
10. 10. 8000. 250. 3 
END 

In spite of the fact that the MATE command must be 
included in this part of the data preparation, its 
discussion is left to the section where the element modules 
are treated. 

After completing the mesh data input, we are ready to 
initiate the problem solution. 

FEAP has modules for variable algorithm capabilities and 
which, if necessary, can be modified or expanded. The basic 
aspect of the variable algorithm program is a macro 
instruction language which can be used to construct modules 
for specific algorithms as needed. This language is 
interpreted by subroutine PMACR and the complete list of 
macro instructions available in this version of FEAP is 
given in Appendix A.3. 

Here, as an example, we will discuss in detail the 
algorithm for the solution of linear steady state problems 
using the original explanation given in Ref. 4. 

To use the macro programming commands, the user only 
needs to learn the mnemonics of the language. If one wishes 
to form the global stiffness matrix the program instruction 
TANG is used (TANG is the mnemonic for a symmetric tangent 


stiffness matrix and for nonlinear elements would form and 
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assemble into the global stiffness the element tangent 
stiffness computed about the current displacement state; for 
linear elements this is just the linear stiffness matrix). 
For a problem with an unsymmetric tangent stiffness the 
macro command UTAN is used. 

If one wishes to form the right hand side of the 
equations modified for specified displacements one uses the 
program instruction FORM. The resulting equations are 
solved using the instruction SOLV. 

Printed output can be obtained using the macro command 
DI SP. 

The above instructions are sufficient to solve linear 
steady state problems, that is, the macro instructions 
TANG ( or UTAN ) 

FORM 

SOLV 

DISP 

are precisely the required instructions to solve any linear 
steady state problem. 

If the same block of instructions is repeated twice or 
more, considerable effort is wasted in preparing the macro 
instruction data, To rectify this, looping commands are 


introduced as the instruction pair 


LOOP n 
NEXT 
which indicates that looping over all instructions between 
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LOOP and NEXT will occur n times. 

Many other classes of problems can be solved using the 
macro instruction list given in Appendix A. With the above 
short explanation and the discussion presented in section IV 
for Heat Transfer problems, the reader will be able to 
construct his own algorithms using the macro programming 


language. 


B. FINITE ELEMENT SOLUTION MODULES 

Most of FEAP is common to a wide range of problems, but 
some of the computations are linked to the type of problem 
to be solved. It is clear that if one wishes to solve an 
Elasticity problem, the calculation of the stiffness matrix 
will be different from that used for a Heat Transfer 
problem, It is also clear that differences also arise if one 
is using triangular elements instead of quadratic 
isoparametric elements as a way of describing the continuum, 

The program was designed such that all computations 
associated with any type of element are contained in an 
element subroutine called ELMTnn where nn is between 01 and 
05 in this version of FEAP. Each element type to be used is 
specified as part of the material property data, following 
macro control MATE, 

Since the objective of this work is the solution of Heat 
Transfer problems, the element modules discussed here are 
the two and three dimensional heat tranfer modules named 


ELMTO2 and ELMTO3. 
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1. Two Dimensional Heat Transfer Element 


Subroutine ELMT02 is accessed using the number two 
in column ten of each material number card. 

According to the value of the variable !SW, defined 
in subroutine PMACR, a specific function is performed. If 
ISw=1 it reads and prints the material property data and 
line boundary conditions. The stiffnes and mass matrices are 
computed if I!SW=3 and !SW=5 respectively. The load vector 
due to internal forces and boundary conditions is calculated 
when | SW=6. 

The material properties may be input as constants or 
as a table of temperature and property values. If any other 
way of defining the property values is required by a 
specific problem, subroutines PKX , PKY , PROC or PO which 
calculate respectively the conductivity in x and y 
directions, the thermal capacity and the heat generation per 
unit volume, can be easily modified or substituted. 

The shape functions used by this heat tranfer module 
are calculated by subroutines SHAPE and SHAP2. These 
routines are capable of constructing shape functions for a 
three-node triangle, a four-node linear quadrilateral, an 
eight-node quadratic serendipity quadrilateral, a nine-node 
quadratic Lagrangian quadrilateral or any combination 
in-between. 

A four- to nine-node two dimensional element is 
shown in Figure 3. As shown by Bathe and Wilson in Ref. 5 


the temperature within the element is expressed at any time 
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in the local coordinate system s, t in terms of the nodal 
temperatures h by 


T(s,t) » x N.(s,t) T. 
i 1 


where 

N, = (1-s)(1-t)/4 - (NetNg)/2 - Ng/4 
N, = (L#s)(1-t)/& - (Ng+Ng)/2 - Ng/4 
№, = (1+5) (1+0)/4 - (Мо+М)/2 - М/% 
Ny = (l=s)(1+t)/4 - (Ny+Ng)/2 - Ny/k 
Nç = (1-s2)(1-t)/2 

E (1+s)(1-t7)/2 

N, = (1-s7)(1+t)/2 

Ng = (1-s)(1-t2)72 

N. 3 (1-52)(1-t2) 


If any of the nodes from five to nine are omitted 
the corresponding value of N is zero. 

The connectivity table must follow the numbering 
convention shown in Figure 3, otherwise wrong results will 
be obtained. 

The maximum number of nodes per element will depend 
on the type of elements to be used and may be from three to 
nine. The eight-node serendipity shape functions occur if 
the ninth node number is omi tted. The four-node 
quadrilateral shape functions are computed if only the first 
four nodal connections are non zero and the three-node 
triangle if the first three nodal connections are non zero. 
Finally, if a mid-side nodal connection Is omitted the edge 


is linear. 
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The line boundary conditions, closely related to the 
type of element used, are defined in the cards following the 
macro MATE. The line subjected to a specified boundary 
condition is identified using the local coordinates. A line 
is numbered 1or 2 if it is perpendicular to the axis s or 
t, respectively. Those numbers are positive or negative 
according to which direction of the axis it is 
perpendicular. Thus, as an example, the line defined by the 
nodes 2, 6 an 3 Is numbered 1 while the line defined by 
nodes 1, 5 and 2 is numbered -2. 

The boundary conditions are treated in subroutine 
BCOND2 and four types are considered: 

a) specified flux 

b) convection with constant heat tranfer coefficient 

c) convection with temperature dependent heat 
transfer coefficient 

d) radiation 

If more complex boundary conditions are required 
subroutine BCOND2 may be easily modified. 

The numerical integration necessary to evaluate the 
matrices and vectors in equation (6) are performed using the 
Gauss quadrature formula; the abscissas and weight 
coefficients are calculated in subroutine PGAUSS. The user 
may choose the number of points per direction used in the 
integration from one to six, but the default value is four 
points . 


To illustrate the above, we assume that our sector 
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used as an example Is part of the cross section of a hollow 
cylinder made of a material with conductivity in x and y 
directions of 10, specific heat 250 and density 8000. If the 
outside surface is subjected to convection to an ambient 
atmosphere at 20 degrees with a constant heat tranfer 


coefficient of 200, the cards following the macro card MATE 


are 

MATE 
1 2 

10. 10. 8000. 250% 3 
6 2 1 200. 20. 


it must be noted that a plane analysis is specified 
and the integration will be performed using three points per 
direction. The reader must also note that this version of 
the problem is slightly different from the one given when 
the mesh data preparation was discussed. In the previous 
problem the line now subjected to a convection boundary 
condition was at a specified constant temperature of 20. If 
this second problem is to be solved the cards corresponding 
to nodes 31-33 in the BOUN and FORC sets must be omitted. 

Following is a list of subroutines included in the 
two dimensional heat tranfer module: 

ELMTO2 : main routine; forms the necessary arrays 
for the solution. 

SHAPE : calculates the shape functions for the 
triangle and the linear quadrilateral elements. 


SHAP2 : adds the quadratic terms and the center node 


22 





to the shape functions, 

PGAUSS : determine the abscissae and weight 
coefficients for the numerical integration. 

BCOND2 : adds to the arrays calculated in ELMTO2 the 
contribuition from the specified boundary conditions. 

PKX : determine the temperature dependent 
conductivity in the x direction. 

PKY : determine the temperature dependent 
conductivity in the y direction. 

PROC : determine the temperature dependent thermal 
capacity. 

PQ : determine the temperature dependent heat 
generation per unit volume. 

CONV : determine the temperature dependent heat 
transfer coefficient. 

TABLE : called by РКХ, РКҮ, PO, PROC and CONV; 
calculates the correspondent coefficient to a gi ven 
temperature by linear interpolation between two consecutive 
entries in a table. 

JACBB2 : calculates the jacobian determinant when 
an integration over a line is necessary. 
2. Three Dimensional Heat Transfer Element 

This module is called ELMTO3 and is accessed using 
the number three in column ten of each material number card. 

This element is the generalization of ELMTO2 for 
three dimensional space and little remains to be said about 


its use. 
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The material properties are read in a similar way as 
the two dimensional case, but the conductivity in the z 
direction must be added. This module uses the same 
subroutines as ELMTO2 to read and calculate the temperature 
dependent properties. 

Subroutines SHAP3D and SHAP3 may construct 
interpolation functions for any combination between a 
eight-node linear brick and a 21-node Lagrangian brick. 

An eight- to twenty one-node three dimensional solid 
element is shown in Figure 4. The natural coordinates (r, s 
and t) of the eight corner nodes are (*1, *1, +1) , of the 
twelve mid-edge nodes are (0, +1, +1), (+1, 0, +1) апа (+1, 
+1, 0) and of the center node are (0, 0, 0). 

The temperature within the element is defined in 
terms of the nodal temperatures T, at any time by [Ref. 5] 

Т(г, 5,1) = У Ni(r,s,t)Ti 

If we define 

GCB, 81) - 5 (148.8) ‚ for 8.2 +1 

= 1-8 for 8,- 0 
and 
g: = GCr,r.2G8Cs,s; )GCt, ti) 


1 


where ri^ S; and t; are the natural coordinates of the 


element nodal points, the interpolation functions N; are 


Ne Ny NN 7212 = Nay Lo 
№, ΞΕ» - (Να “ΝΑ ΗΝ α)/2 - Νρη/8 
Ν = g3 7 (N¿p+N7]+N]92/2 - N21/8 
N4 BUE. = (N44 *N15*N59) /2 - N51/8 
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Ns = 55 eS (Ni ;*Ni G*N4 72/2 в №, 1/8 
Nç = £g - (Nig5*N,4*N4g)/2 - Νρ/8 


Ент о 5578 

Ng = gg ~ SNyctNy otha) /2 - N. /8 

N, = g; = Nz7/4 for j=9,...,20 
a Pol 


If any of the nodes from nine to twenty one are 
omitted the corresponding value of g is zero. 

The numbering convention used in the definition of 
the shape functions and shown in Figure h must be followed 
by the user when the connectivity table is built and the 
surface boundary conditions are specified. The surfaces are 
numbered +1, +2, +3 as they are perpendicular to the 
positive or negative directions of T axis r, s, t, 
respectively. 

As the numerical integration routine is the same as 
for ELMTO2, the user may choose from one to six points per 
direction, the default value being four points. 

Following is a list of subroutines included in the 
three dimensional heat transfer module: 

ELMTO3 : main routine; forms the matrices necessary 
for the solution. 

SHAP3D : calculates the shape functions for the 
eight node linear brick. 

SHAP3 : adds the quadratic terms to the shape 
functions and the center node for the Lagrangian brick. 


BCOND3 : adds to the matrices the contribuition from 
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the specified boundary conditions. 
PKZ : determine the conductivity in the z direction 
JACBB3 : calculates the jacobian determinant when an 


integration over a surface is necessary. 


C, TIME INTEGRATION MODULE 

The first order ordinary differential equation solver is 
accessed when the macro command ODEI is utilized and uses 
the Zienkiewicz two- and three-time level schemes. The 
details of these algorithms are treated in section IV. 

This module was designed for the solution of first order 
equations but can, with little programming effort, be 
expanded in order to include higher order algorithms and 
solve higher order differential equations. 

In order to economi ze computer memory space at some cost 
of execution time, this module only uses the space 
previously reserved for the mesh construction. As the mesh 
data is not needed during the execution of a time step 
integration, all the corresponding arrays are kept in file 9 
during that period. Once the next displacement vector is 
calculated and the current displacement vector is saved in 
file 10, all the information in file 9 is retrieved and the 
mesh data restored. 

File 10 is used as a working space during the time 
integration and the current displacement vector is saved in 
it for possible future use. This vector is retrieved if the 
algorithm used in the next time step integration is a three 


point scheme. 


26 





Besides the time step size change using the appropriate 
macro instructions (macro TOL), an optional automatic time 
step adjustment was incorporated in the program. The norm of 
the difference between the displacements vectors at two 
consecutive times is computed at each step. If the norm is 


less than a predetermined value AT the time step size 


max ^4 
is doubled before going to the next step, whereas if the 


norm is greater than a prespecified value AT the time 


min’ 
step size is halved and the calculation for that time step 
is repeated until the norm is acceptable. The magnitudes of 
the maximum and minimum values of the norm are problem 
dependent and must be choosen by the user. If they are not 
specified no time step adjustment will be performed. The 
recalculation of the displacement vector is allways done 
using the two point scheme, even when a three point scheme 
was utilized for the first calculation. 

The macro command ODEl used to access this module 
must be followed by a second macro command in order to 
determine the function to be performed. The secondary macro 
instructions are INIT, LINE or QUAD. 

The couple ODE1 INIT reads the integration constants 
theta , beta and gamma, the parameters for the automatic 
time step adjustment and the initial displacement vector. 
This data must follow the macro program (see user's 
Instructions). No time integration is performed by this 
instruction. 


The couple ODE1 LINE performs the two point scheme 
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and the current displacement vector is substituted by the 
new lcılated displacement vector. 

The couple ODEl QUAD has a function similar to ODEl 
LINE but uses the three point scheme. 

Once the execution of an ODEl macro command is 
complete, the mesh data is restored and the displacement at 
time t replaced by the displacement at time t+ At, when the 
secondary macro command is LINE or QUAD. If INIT is used, 
the displacement vector becomes the given initial 
displacement vector. 

The subroutines included in this module are: 

PODE1 : control routine; reserves working space and 
calls the appropriate subroutines. 

INIT : reads the constants theta, beta and gamma and 
the maximum and minimum values for the norm used in the time 
step adjustment. 

PLINE : performs the two point integration scheme. 

PQUAD : performs the three point integration scheme. 


TERM : performs the automatic time step adjustment. 
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IV. APPLICATION OF FEAP TO HEAT TRANFER PROBTEMS 


The program FEAP requires from the user the knowledge of 
the algorithm to be used in the solution of the problem to 
be treated. 

In this section the algorithms used in the solution of 


heat conduction problems are discussed. 


A. STEADY STATE PROBLEMS 
1. Algorithms for Linear Problems with Macro Program 


In this case equation (8) becomes 


[K] (T) + {F} = {0} (12) 
The algorithm used to solve (12) is described by 

[K] fv} = -{F} -[K] {T°} (13) 

(T) = (T°) + (v) (14) 


where {T°} are the specified nodal temperatures. 

The macro instruction TANG builds the left hand side 
of (13) while FORM forms the vector in the right hand side. 
The instruction SOLV solves the system (13) and performs the 
step defined by (15). 

Consequently the simplest macro program used to 
solve this problem is 
TANG | 
FORM 
SOLV 
DISP 


29 





2. Algorithms for Nonlinear Problems with Macro 
Program 
Equation (12) still applies to this case but now [K] 
and {F} may be dependent on the nodal temperatures {T}. The 
well known Newton-Raphson iteration may be used and the 
algorithm is described as 
[K(T*)] (vi) = -{F(T*)} - [K(T?)] (тър үк (15) 
S E o (16) 
where (T are the nodal temperatures at the ith iteration. 
The algorithm defined by (15) and (16) is the one used for 
the linear problem repeated several times. Then the macro 
program may be 
LOOP n 
TANG 
FORM 
SOLV 
DISP 
NEXT 
DISP 
where n are the user's guess of the number of iterations 


necessary to obtain equilibrium. However, the program has an 


internal check on the value of the vector norm ||β7|| , 
where 
m2 1/2 
IRI] = (5 [А 19 
k=1 
Whenever 
IIRTI]< TOL x max | IR || (j=1,...,1) 


J = 
where TOL Is a predefined tolerance ( 17° is the default 
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value), the iteration ceases and a skip to the macro command 
immediately following the first NEXT occurs. Usually one 
begins with zero as the initial guess of the displacement 


vector; however, any other vector may be used. 


B. TIME DEPENDENT PROBLEMS 
1. Two and Three Point Recurrence Schemes for First 

Order Equation 

The set of differential equations 

[K] (T) + [C](T) + (F) = (O) (8) 
may be solved using one of the many recursive schemes 
described in the literature. From them, the two and three 
time level schemes presented by Zienkiewicz [Ref.4] offer a 
wide range of choices for the solution of linear and 
nonlinear problems. 

The recurrence relation for the two-time level 
scheme can be written as 

( 1С] + өГК]) (Т, 1} + ELO + (1-9)[K] (T ) + 

{F} = {0} Е 0 < 0 < 1 (17) 
where the subscripts n and n+1 denote evaluation at time t 
and t+At. The vector (F) is defined as 

ЕЕЕ r) (18) 

The choice of the parameter @ defines the particular 
scheme to be used and the reader will recognize a well known 
series of finite difference formulas with a modification of 
using a weighted loading term {F}. 

Consider the system of decoupled equations in terms 


of the modal participation variables T; 


D 





СЕЕ = 0 (19) 


for free response, i.e., ia, expresston (17) becomes 


(+ €; * k19) (Ti)444 * + c, + k,(1-0)) (T,) =0 (20) 
Substituting the relation 
(o my = MTI (21) 
in (29), the resulting characteristic equation of the 
recurrent scheme solved for À gives 
À = [1-k, (1-0)At/c,]/[1+k,0 At/c;] 
The scheme will be unconditionally stable if 
a] <1 
for any value of p; where 
р; = (К. /с;) At (22) 
This condition is satisfied for 
© > 1/2 
On the other hand if 
0«0« 1/2 
stability is conditional requiring 
p; < 2/(1-20) 

Figure 5 shows how A varies with Pi for the schemes 
discussed later and how it compares with the exact value of 
\ = ехр(-р;) 

The schemes considered are 

a) Crank-Nicolson with 0=1/2 

b) 0=2/3 proposed by Zienkiewicz [Ref. 4] 
c) 025/h proposed by the author 

d) 020.8783 proposed by Liniger in Ref. 6, 


The Zienkiewicz three-time level scheme is defined 
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as 
(Y[C] *8 At[KD O44). * ((1-2y)[C] + (1/2-2B+Y)At[K]) (T), )+ 
(-(1-Y)[C] + (1/2+8-y)At[K])tT, _,) NEO (23) 
where the subscripts n, n+l and n-1 denote evaluation at 
time t, t+At and t-At respectively. 

The vector [ F) is defined as 
ΗΠ. : ΠΥ2-26»γ}15.} - (1/2«B-y)tF, .) (24) 
and the parameters y and ß define the particular scheme to 
be used. 

Considering again the system of decoupled equations 
(19) and using the relations (21) and (22) one finds that 
the characteristic equation for the three-time level scheme 
is 
(y+8p,) AS + [(1-2y) + (1/2-28+v)p,]2 + 

[-(1-Y) + (1/2+B-Y)p;] = O (25) 

Writing 

g = [1+(1/2+y)p,]/[y+6p;] 
h = [-1+(1/2-y)p,]/[y+8p; ] 


the roots of (25) may be written as 

му о = (2-g)/2 + [(2-g)2-4(1-n)]1/2/2 

These roots will be complex if the quantity under 
the square root is negative. Then the modulus of À is 

la] = (1-1)1/2 

The scheme is stable if |A|<1 for any р; апа Wood in 
Ref. 7 shows that this condition is satisfied if 


Y>1/2 and В>ү/2 
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Equation (25) has two roots: Az + the principal root 

which is the aproximation to the exact 
\ = ехр(-р;) 
and the spurlous root hos For relative stability one must 
have 
ТАБИ 

Also а negative root or complex roots can produce 
oscillation. 
and A, and the 


1 2 
modulus |, | for the:complex roots are plotted against p, 


In Figure 6a-f the real roots A 


for the schemes: 

a) B=1/3, y=1/2 proposed by Lees in Ref, 8 

b) B=3/4, y=1 proposed by Hogge in Ref.9 

c) 830.656, y=1.2184 proposed by Wood [Ref. 7] 

а) 8=4/5, y=3/2 proposed by Zienkiewicz [Ref. 4] 

e) ß=9/10, y=3/2 proposed by the author 

f) Fully implicit algorithm, 8=1, y=3/2 [Ref. 7]. 

From Figure 6a one may predict a very strong 
oscillatory behaviour for the Lees algorithm a). 

In order to study the various schemes, the physical 
problem represented by the nondimensionalized linear heat 


conduction equation 


ET 
st Sx? 


was solved in a bar of length 4 and width 1 subjected to the 


boundary conditions 


T=1 at x=0 
eT. At x=4 
6X 
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and initial condition 
i T=0 at t=0 

i.e., consider a step change in the surface temperature. 
This problem was also studied by Wood and Lewis in Ref. 10. 

The problem was discretized in space using ten equal 
length two dimenstonal quadratic elements through the length 
of the bar. This type of elements was chosen because it was 
reported by Wood and Lewis as giving the worst numerical 
results. 

in order to avoid the discontinulty in the loading 
term caused by the step change of the surface temperature at 
the Initial time, the problem was first solved starting from 
time t=At, where At Is the time step size. The value of the 
temperature distribution at t=At is provided by the exact 
analytical solution. This procedure also provides the 
necessary two starting vectors for the three-time level 
schemes. 

The temperature of the nodes at x=l is used as 
reference value to compare the various schemes. 

in Figure 7a-d the results obtained from the 
various two-time level schemes with At=2 for the temperature 
at x=l are compared with the analytical solution. For this 
Ideal starting conditions the Crank-Nicolson, 9=1/2, proves 
to be the most accurate algorithm. When the step size was 
Increased to 10, al] the schemes performed well as shown in 
Figure 3a-d. 


The corresponding results for the three-time level 
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schemes are presented in Figure 9a-f and Figure 10a-f for 
At=2 and At=10, respectively. 

Strong oscilations were observed with the Lees 
algorithm a). 

For At=2 the schemes c), d) and e) performed well. 
For the larger step size of 10 the safest schemes are e) and 
D. 

The step change In the loading term has been 
reported [Ref. 10] as producing relative instability when 
some of the schemes are used, particularly with the 
Crank-Nicolson algorithm. The noise may be supressed 
reducing the magnitude of the time step but this is 
impraticable In most problems. Other techniques may be used 
to reduce the amplitude of the oscillations as those 
discussed by Wood and Lewis [Ref. 10] and Gresho and Lee in 
Ref. 11. 

The step change In the loading term in the problem 
in study Is due to the variation of the surface temperature 
at the initial time t=0. That discontinuity In the force 
term is illustrated in Figure lla. When the time dimension 
is discretized it has been common practice to calculate the 
numerical solution at the nodes t=0, t=At, etc. When a 
three-time level scheme is used the node at t=- At is 
considered assuming all conditions steady before t=0. The 
force values used by this method are indicated by the 
circles in Figure lla, This procedure transfers to the 


numerical solution the uncertainty In the value of the force 
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term at time t=0 introduced by the mathematical 
discontinuity assumed In the analytical solution. This 
problem is avoided if one uses the procedure illustrated in 
Figure 110. The nodes considered for the discretization are 
those at t=-At/2, t=At/2, t=3At/2, etc. The discontinuity at 
t=0 is smoothed by the interpolation of the loading term 
inherent to the particular scheme being used. The values of 
the force term for every node are well defined and no 
uncertainty is associated with any of them. For the 
three-time level schemes the node at t=-3At/2 is used 
assuming all conditions steady before t=0. 

The problem studied was solved using both starting 
procedures. The results obtained with the first procedure, 
start at t=0, using the two- and three-time level schemes 
are shown in Figures 12a-d and ]13a-f, respectively. The 
second procedure, start at t=- At/2, gives the results shown 
in Figures lha-d and 15a-f for the two- and three-time level 
schemes, respectively. 

One may conclude that, in this simple problem, when 
the procedure Illustrated in Figure 11b is followed, the 
step change in the loading term do not deteriorate the 
accuracy of the numerical results obtained with the schemes 
tested, with the exception of the Crank-Nicolson algorithm. 
This scheme shows an oscillatory behaviour not observed when 
the step change In the force term is not present. 

The starting procedure of Figure llb also proves to 


be an efficient way of starting the time integration with 
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the three-time level schemes, 


2. Algorithms for Time Dependent Problems with Macro 
Program 


A simplified form of the two-time level scheme is 
available in FEAP and is described by 
AC] + OK l) Tae} + CALC] + (1+6) [K,]){T,}+ 
At{F} = {0} (26) 
where the subscripts п апа п+1 denote evaluation at time t 
and t+At, respectively, 

The three-time level scheme programmed in FEAP is 
described by the expression 
(y[C,] +8 At(K,]){T,,,} + ((1-2y)[C, ] + (1/2-28+y)at[K,]) {7} 
* CO- IC, » (1/2+8-y)At [K ]) {T 1? 

+ {F} = {0} (27) 
where the subscripts n, n*l and n-1 denote evaluation at 
times t,At* t and t-At respectively. 

It must be noted that in (26) and (27) the vector 
USE an aproximation of the Interpolated force vector {F}. 

If the force vector Ís strongly dependent on time, 
this aproximation may produce wrong results. 

In problems where the stiffness matrix is constant 
and the force term results from the specified boundary 
displacements or where the specified boundary forces are 
only time dependent, the interpolation of the force term may 
be easily acomplished since the force/displacement boudary 
values can be changed at any time using the macro commands 


MESH or PROP. This procedure was found partcularly useful 
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when a step change in boundary displacements or forces is 
applied at the initial time, 

The first order ordinary differential equation 
solver module is acessed using the macro command ΟΡΕΙ 
followed by one of the following macro commands: 

INIT to specify the initial displacement vector 

and the integration constants 

LINE to perform the two-time level algorithm 

QUAD to perform the three-time level algorithm 

Since the matrices [K] and [C] and the vector {F} 
must be formed before the time integration, the macro 
commands TANG, CMAS or LMAS and FORM are closely associated 
with the use of ODEl. It should be noted that the macro FORM 
forms the vector 

{RCT _) } = EL - [Ka] {Τ} 
which is allways dependent on the current displacement. 
Therefore the macro command FORM must be used every time 
step and preceed ODEI. 

To solve a fully nonlinear problem the following 


macro program may be used: 


DT At 
ODEl INIT 
LOOP n 
TANG 


CMAS ( or LMAS ) 
FORM 
ОрЕ1 LINE ( or QUAD ) 
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TIME 
DISP 
NEXT 

The macro command program must be followed by the 
data cards necessary to specify the integration constants 
0, Y and B and the initial displacement vector. 

If the matrices [K] and/or [C] are not temperature 
or time dependent then the instructions TANG and/or CMAS or 
LMAS must be placed outside the loop. 

For the simplest case of a linear problem, the macro 
program may be 
DT At 
ODE] INIT 
TANG 
CMAS ( or LMAS ) 

LOOP n 

FORM 

ODEl LINE ( or QUAD ) 
TIME 

DISP 

NEXT 
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V. NUMERICAL EXAMPLES 





A. RADIAL TEMPERATURE DISTRIBUTION IN A HOLLOW 

CYLINDER 

Consider the hollow cylinder of infinite length already 
discussed in section III. Consider the case in which the 
inside wall is maintained at a constant temperature and heat 
is lost by convection through the outer wall. This problem 
was considered for solution by Lew [Ref. 3]. 


The geometrical and thermophysical characteristics are: 


Inside radius (rig) 0.1 m 
Outside radius (rout? 9.3 m 
Thermal conductivity (k) 10 W/m °C 
Outside ambient temperature 20726 
Heat transfer coefficient 200 W/m**C 
Inside wall temperature (Ti) 820 °C 


The finite element model of this problem was completly 
discussed before and is shown in Figure 1. The complete 
computer output is given in Appendix B. It must be noted 
that the element arc width used is arbitrary since the 
problem is symmetric and no heat is conducted perpendicular 
to the radius. 

_ The comparision between the analytical and FEAP 
solutions is shown in Figure 15. In this case the values 


provided by the finite element model used are very accurate 


and no further mesh refinement is necessary. 
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B. NONLINEAR STEADY STATE HEAT CONDUCTION IN A SLAB 

WITH TEMPERATURE DEPENDENT CONDUCTIVITY 

A liquid is boiled by a flat electric heater plate of 
thickness 2L. The internal energy Q generated electrically 
may be assumed to be uniform. The boiling temperature of the 
liquid, corresponding to a specific pressure, is Tye 

In the finite element analysis for the temperature 
distribuition in the heater, ten equal length quadratic 
elements were used to represent the unit cross section 
through the half plate thickness L, since the problem is 
symmetric. The extremity of the resulting slab is assumed at 
the temperature T. while the other is insulated. 

The thermal conductivity is assumed temperature 
dependent according to the expression 

k = a(1+bT) 

where a and b are constants. 

Assuming the following values for the parameters іп a 


consistent system of units 


T = 0,0 
W 

а = 0,5 
L = 1.0 
Q = 1.0 


the problem was solved for b equal 0.0, 1.0 and 2.0. 
The numerical results are compared with the analytical 
solutions in Figure 16 and one may conclude that FEAP 


provides exact solutions in this problem. 
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C. TRANSIENT TEMPERATURE DISTRIBUITION IN A HOLLOW 

CYLINDER 

The same hollow cylinder of infinite length treated 
before is considered again but the case solved is the one 
with both outside and inside walls maintained at constant 
temperatures. Thermal diffusivity, a=k/ pc is specified as 
0.00005 ee Initially the cylinder is at an ambient 
temperature Τρ of 20°C: suddenly the temperature of the 
inside wall is raised to 820°C while the outside wall is 
maintained at 20°C. All other conditions are the same as for 
the steady state problem. 

In order to obtain the unsteady solution the number of 
elements in the finite element model was doubled, thus 
twelve radial quadratic elements were used. The element arc 
width chosen is the same 6° as before. 

The two time-level scheme with 9=2/3 was used for the 
time integration. In order to compare the results the 
dimensionless temperature T*, radius r* and Fourier number 


Fo were used. They are defined as 


х — - - 
T (T TT, To) 
х = 
n г 
апа 
u 2 
Fo = ατ/τ{ῃ 


A time step size At of 2 seconds was chosen initially. 
After 10 sec, At was increased to 10 sec and to 50 sec after 
200sec. After 1000 sec,At was increased to 100 sec. 


Figure 17 shows the evolution of the temperature of the 
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center interior points of the cylinder and Figure 13 
compares the analytical values with the FEAP solution for 
the radial temperature distribuition for two different 


times. 


D. SLAB WITH RADIATION-CONVECTION BOUNDARY CONDITIONS 

Consider a plane slab of thickness L, with constant 
thermophysical properties (ks p c=1), subjected to 
simultaneous radiation and convection at x=0 and insulated 
at x=L (0<x<L). The slab is initially at a uniform 
temperature Τρ and it was decided to consider zero ambient 
temperatures for convection and radiation. 

lt was also decided to consider the Biot number 

Bi = h L/K = 1 
where η. is the convection heat tranfer coefficient and 

R = eoToL/k = 4 
where e is the emissivity and c the Stefan-Boltzman 
constant. 

The same mesh as in Example B was employed for the 
finite element solution. The time integration was performed 
by the three-time level scheme with y =1/2 and f=1/3. A 
variable time step size was chosen. A value of At=0.001 was 
used initially. After t=0.02, At was increased to 0.0025. 
After t=0,2, At was increased to 0.01. After t=0.4, At was 
increased to 0.025. After t=1.0, At was increased to 0.05. 

The history of the ratios T/T and Τι / Τρ are plotted in 
Figure 19 against the Fourier number 


BT 
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- is the surface tmperature at x=0 and Um is the surface 
temperature at x=L. 

The finite element results are compared with those given 
by Haji-Sheick and Sparrow in Ref. 12 and obtained froma 


Monte Carlo method. 


45 





VI. CONCLUSIONS. AND RECOMMENDATIONS 


The computer program in this thesis provides an accurate 
and reliable means for solving a variety of Heat Transfer 
problems. The use of this program and efforts to increase 
its versatility are highly encouraged. 

The capabilities of the program, while qui te 
significant, can still be improved. The present version is 
designed for an "in-core" solution technique, which 
restricts the problem size to within the computer core-size. 
The capacity to handle large problems may be increased 


"out of core" technique for solving 


through the use of some 
the system of equations and even for building the mesh data. 
In that case the size of the problems treated would be 
restricted only by the availability of external storage 
devices. 

In order to check the mesh input and to easily interpret 
the output, a graphics module must be included in the 
program. 

Although the time integration module seems to provide 
reliable results in linear and some nonlinear problems, it 


should be subjected to further research in order to test it 


with different types of nonlinear analvsis. 
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Figure 1. Mesh for Hollow Cylinder Problem 





Figure 2. Eight-Node Serendipity Quadrilateral 
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Figure 3. Local Node Number Sequence for a Quadratic 
Lagrangian Element 
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Figure 4. 


Local Node Number Sequence for a 21-Node 
Three Dimensional Element 
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